1. Field of the Invention
Disclosed is a method for calculating static corrections for seismic data by minimizing an objective function. The surface consistent source and receiver statics that maximize the coherence of reflections are determined by simulating the crystallization process of a system of fictitious classical particles whose interaction potential is given by the objective function.
2. Discussion of Related Art
Seismic exploration provides a means for imaging in great detail, the structural configuration of subsurface earth layers. As is well known, a plurality of several tens or several hundred seismic receivers are emplaced along a desired line of survey at spaced-apart intervals as a receiver spread. An acoustic source generates a wavefield at or near the surface of the earth. The wavefield propagates downwardly and impinges upon subsurface strata whence it returns to the surface where the reflected wavefield is detected by the seismic receivers which convert the seismic signals to varying electrical signals. The detected seismic signals are transmitted to a multichannel recording device over an electrical, optical or ethereal data communications link. The received signals are recorded on an archival storage medium of any desired type as a data set comprised of a plurality of raw, usually discretized, time scale recordings, one recording for each receiver or receiver group.
After suitable data processing, the members of the raw data set are converted to visual form to be displayed on a cross section as a series of side-by-side variable-amplitude or variable area time-scale oscillographic traces. The section depicts a two-dimensional cross section of the earth having the dimensions, vertically, of two-way travel time in seconds and spatially, along the horizontal axis, as a function of the receiver spacing as measured in desired units such as meters. Typical dimensions for a cross section created from a single data set might display reflection events embraced within a vertical travel time window of four to eight seconds or more, along a horizontal spread corresponding to fifty to one hundred receivers or receiver groups spaced, perhaps, 100 meters apart.
Each firing of the acoustic source is followed by a listening period of a desired length during which the reflected wave fields are recorded. The source and the receiver spread are moved along the line of survey by a multiple of the receiver spacing and the source is fired again to record another data set. That process is repeated until the entire length of the line or area of survey has been occupied.
Wind, vehicular and pedestrian traffic, microseisms and at sea, ship noise and environmental noise due to marine life, may seriously interfere with and obscure the desired reflection data. Furthermore, solid friction in the earth as well as wavefield scattering attenuate the reflected signals. For that reason, a multiplicity of traces from different shots resulting from the interaction of a plurality of cooperating sources and receivers, but having the same subsurface incident point, are gathered and combined together by use of various algorithms well known to the geophysical profession. Multiple coverage tends to destructively attenuate random noise but enhance coherent reflection data. One such method is referred to as a Common Mid Point stacking which will now be illustrated with reference to FIGS. 1-3.
In FIG. 1, acoustic source $3 located on surface 10 insonifies a subsurface earth layer such as 14 whence the wavefield is reflected therefrom along ray paths 16, 18, 20 to be detected by receivers R1, R2, R3. The reflected wavefield as seen by each receiver is recorded either in analog or digital format on a separate recording channel. The wavefield as recorded at each receiver station may be presented visually as plurality of oscillographic time scale traces, one trace per channel. The resulting exemplary set of three records represent a common source gather.
In FIG. 2, sources S1-S3 respectively insonify receiver R3, after wavefield reflection from interface 14 along ray paths 20, 22, 24. The resulting recordings represent a common receiver gather. Of course, sources S1-S3 also insonify receivers R1 and R2 but those ray paths have not been shown to avoid complicating the drawings.
Observe that the subsurface incident points are different for each source-receiver combinations of FIGS. 1 and 2. Assuming that interface 14 is flat, the three records from either the common source gather or the common receiver gather or all six records from both gathers could be combined or stacked together after correction for angularity (also called Normal Moveout or NMO), into a single, zero-offset relatively noise-free trace. In effect, that stacked trace is a raw gather created from the interaction between a plurality of cooperating seismic sources and receivers.
In actual field conditions interface 14 is not necessarily flat. Therefore a preferred process employs common mid-point stacking as now shown in FIG. 3.
In FIG. 3, ray paths are shown for source positions S1-S3 with respect to receiver positions R1-R3. The ray paths shown all converge at a common mid point 15. The surface separation between the respective sources and correspondingly numbered receivers is termed the offset. As can be seen readily, application of NMO to the respective ray paths would collapse them to zero-offset as represented on the surface by S0/R0 to form two-way travel path 26 impinging on common mid point 15. In this case of zero dip, mid point 15 also forms a common depth point as well. Since all of the ray paths are incident on the same subsurface point, the traces can be stacked validly to enhance the signal-to-noise ratio.
Near-surface lateral velocity variations and surface elevation changes create travel-time variations that may be approximated by surface-consistent static time shifts. In FIG. 3 there are shown two subsurface earth layers, separated by dashed line 12, characterized respectively by acoustic propagation velocities of V.sub.0 =1500 meters per second (m/s) and V.sub.1 =3100 m/s. P.sub.si and P.sub.rj indicate those portions of the near-vertical ray paths that traverse the variable-thickness upper low velocity layer, after ray bending due to Snell's law of refraction at the interface, beneath sources S.sub.i and receivers R.sub.j. Because of a longer combined path length of P.sub.s3 +P.sub.r3 through the low velocity layer, the total travel time for ray path S3-R3, after application of NMO, will be longer than the total travel time along a ray path such as S1-R1, corrected for NMO, where the combined ray path P.sub.s1 +P.sub.r1 through the low velocity layer is shorter. It should be understood that under actual field conditions, even though the low-velocity segments of the total ray path are relatively short, travel-time differences are not negligible due to the very low velocity at or near the surface. The respective time delays s.sub.i and r.sub.j due to a variable-thickness low velocity layer are defined as the surface-consistent statics which must be applied as corrections to the reflection travel times prior to stacking for maximizing inter-trace reflection coherency as will be shown in connection with FIG. 4.
The relative thickness of the low velocity layer beneath the respective source and receiver stations is not necessarily known a priori. Therefore, the static corrections are estimated from the pattern of correlatable reflected events as seen on a set of a limited number of adjacent raw seismic traces. In FIG. 4, by way of example but not by way of limitation, there are shown three raw zero-offset seismic traces 28, 32, 36, each with a correlatable reflected event 30, 34, 38 respectively, corresponding to the two-way wavefield travel times between sources S.sub.i and receivers R.sub.j of FIG. 3. Using the reflection 30 on trace 28 as reference, the reflection pattern shows that reflection 34 is delayed by a static time difference equal to .DELTA.t.sub.1 and reflection 38 of trace 36 is delayed by .DELTA.t.sub.2, the relative delay times being due to the different ray-path lengths through the low velocity zone above interface 12 of FIG. 3 as explained above. The statics are termed surface-consistent because they are due to irregularities of the near-surface low velocity layer.
In the trivial example above, the relative static corrections were determinable by inspection. In the real world, in the presence of many thousands of seismic traces, many of indifferent or poor quality, the seismic traces are processed automatically by computer, by minimizing the root-mean-square error between the fitted travel times and the measured travel times in the least-squares sense. Other, more sophisticated methods employing simulated annealing based on Monte Carlo procedures or other generic algorithms have been proposed to improve the estimates and to minimize expensive computer processing time.
D. H. Rothman, in a paper entitled Nonlinear Inversion, Statistical Mechanics and Residual Statics Estimation, Geophysics, v. 50, n. 12, 1985, reformulates the residual statics problem as an optimization problem. He implemented a standard simulated annealing scheme based on Monte Carlo procedure in order to maximize the stack power which is a measure of the quality of a stack. He chose the random statics as the initial statics of the data. For a given station, a new random static is assigned and the stack power is recomputed. If the new static shift gives a higher stack power, the static shift is accepted for the given station. If the new static shift decreases the stack power, a random number is generated and compared with exp(-.DELTA.E/T) where .DELTA.E is the change in stack power and T is the temperature of the system. If the random number is greater then exp(-.DELTA.E/T), the static is rejected. Otherwise the static shift is accepted as the new static although it decreases the stack power. The same procedure is repeated for all of the stations to complete an iteration. The iterations are repeated with appropriated annealing schedules, i.e., decreasing temperature. As the temperature decreases, exp(-.DELTA.E/T) decreases correspondingly and so does the acceptance rate. As the name suggests, the evolution of the Monte Carlo approach relies solely on luck (random numbers). Of course, a Monte Carlo approach will find the global maximum stack power if we try enough random numbers. But that requires an enormous amount of computer time.
In the same issue of Geophysics, Joshua Ronnen and J. F. Clarebout present a paper Surface Consistent Residual Statics Estimation by Stack Power Maximization. Here, each station is considered one at a time. First, all of the static shifts, except that for the current station, are fixed, then a search is made for that static shift which maximizes the stack power. The procedure is repeated for all of the stations, Unfortunately, that method frequently does not lead to a global maximum but leads to one of many local maxima.
Another drawback to the above methods is that stack power is not a good objective function. that is, a function that is to be maximized or minimized. As Ronnen and Clarebout showed in their paper, the arbitrary linear trend in statics will not change the stack power. As a result, the statics solution which maximizes the stack power, frequently gives a bad stacked section and it is almost impossible to remove a linear trend.
In a later paper, published in Geophysics, v. 51, n.2, 1986, Rothman published a modification of his earlier method. He uses explicit cross correlation of traces instead of simply picking peaks of the cross-correlation functions and transforms the cross-correlation functions to probability distributions and then draws random numbers from the distributions. Estimates of statics are iteratively updated until convergence to the optimal stack is achieved.
The methods cited above are expensive of computer time. There is a need for a more economical deterministic computation method for estimating the residual statics for massive numbers of seismic traces.